Quantum Particle in a Box

A little paradox, traced to a failure of hermiticity

Nicholas Wheeler
10 January 2010
Revised 16 January, in response to some observations by David Griffiths, Darrell Schroeder & Joel Franklin

Introduction

A quantum particle of mass m moves on the interval 0 < x < a. The Hamiltonian opearator is standardly taken to have the form

As a matter of computational (mainly graphic) convenience I set a = 1. The normalized energy-eigenstates are

e[x_,n_]:=(Sqrt[2])Sin[n π x]
Table[(∫_0^1)e[x,m]e[x,n]\[DifferentialD]x,{m,1,4},{n,1,4}]//TableForm

Plot[{e[x,1],e[x,2], e[x,3],e[x,4]},{x,0,1}, Ticks->{{1},{-(Sqrt[2]),(Sqrt[2])}}, PlotRange->{-2,2}]

From

ℰ=(((π^2) (ℏ^2) )/(2 m));

-((ℏ^2)/(2m))D[e[x,1],{x,2}]==ℰ e[x,1]
-((ℏ^2)/(2m))D[e[x,2],{x,2}]==(2^2)ℰ e[x,2]
-((ℏ^2)/(2m))D[e[x,3],{x,2}]==(3^2)ℰ e[x,3]
-((ℏ^2)/(2m))D[e[x,4],{x,2}]==(4^2)ℰ e[x,4]

we see that the (n^th) energy eigenvalue is given by

The paradoxical results that are the point of this discussion emerge when we assume the (momentary initial) state of the system to have the simple form of the following…

Parabolic test state

We declare an interest in the polynomial state

ψ[x_]:=(Sqrt[30]) x(1-x)

which is seen to be normalized:

(∫_0^1)((ψ[x])^2)\[DifferentialD]x

I use a graphic technique to compare the parabolic state to the ground state:

Plot[{ψ[x],e[x,1]},{x,0,1}, Ticks->{{(1/2),1},{(Sqrt[2])}}, GridLines->{{0.5},{(Sqrt[15/8]),(Sqrt[2])}}]

In constructing the figure I have taken note of the fact that ψ[x] is maximal at x = 1/2, where it assumes the value

ψ[(1/2)]==(Sqrt[15/8])

The figure establishes that the ground state is by itself a good approximation to ψ. 

Fourier analysis enables us to improve upon the approximation:

Table[{n,(∫_0^1)ψ[x]e[x,n]\[DifferentialD]x},{n,1,13}]//Simplify//TableForm

That—as these results display—only the odd eigenstates contribute to the development of ψ is made obvious by an elementary symmetry argument. A very simple formula serves to describe the Fourier coefficients of odd order:

Table[{n,(∫_0^1)ψ[x]e[x,n]\[DifferentialD]x-((8 (Sqrt[15]))/((π^3)(n^3)))},{n,1,13,2}]//Simplify//TableForm

We are led thus to define expansion coefficients

c[k_]:=(1/(((2k+1))^3))((8 (Sqrt[15]))/(π^3))

which are found to satisfy the normalization condition

(∑_(k=0)^∞)((c[k])^2)

and manage to do so in excellent approximation even when the series is severely truncated:

Table[{n,N[(∑_(k=0)^n)((c[k])^2),12]},{n,0,10}]//TableForm

With their aid, we construct

(∑_(k=0)^∞)c[k]e[x,2k+1]

That function does not look to be real-valued, and certainly does not look at all like the simple polynomial (Sqrt[30])x(1-x), but seems on the following graphic evidence to check out:

SpectralΨ[x_]:=-(1/(π^3))I (Sqrt[(15/2)]) (E^(-I π x)) (-LerchPhi[(E^(-2 I π x)),3,(1/2)]+(E^(2 I π x)) LerchPhi[(E^(2 I π x)),3,(1/2)])

Plot[{ψ[x],SpectralΨ[x]},{x,0,1},Ticks->{{1},{N[(Sqrt[15/8])]}}, GridLines->{{},{N[(Sqrt[15/8])]}}]

The preceding Fourier series provides already quite a good-looking approximation even when one abandons all but the first few terms:

Plot[{ψ[x],(∑_(k=0)^5)c[k]e[x,2k+1]},{x,0,1},Ticks->{{1},{N[(Sqrt[15/8])]}}, GridLines->{{},{N[(Sqrt[15/8])]}}]

How good? Look to the difference:

Plot[{ψ[x]-(∑_(k=0)^5)c[k]e[x,2k+1]},{x,0,1}]

Note the occurance of Gibbs phenomenon, which becomes progressively more marked (relatively, not absolutely) as the order of the approximation increases:

Plot[{ψ[x]-(∑_(k=0)^20)c[k]e[x,2k+1]},{x,0,1}, PlotRange->{-.0002,.0002}]

Expected value of energy

Let the Hamiltonian

be abbreviated

where evidently

ℋ=((ℏ^2)/(2m));

From

((-ℋ(∂_(x,x))e[x,n])/(e[x,n]))

we recover the information that the (n^th) eigenvalue can be written ℰ (n^2) with

ℰ== (π^2) ℋ

Working directly from

we get

ExpectedEnergy=(∫_0^1)ψ[x](-ℋ D[ψ[x],{x,2}])\[DifferentialD]x

If, on the other hand, we display ψ as a weighted superposition of energy eigenfunctions (i.e., if we work from the spectral representation of ψ) we get

SpectrallyExpectedEnergy=ℰ(∑_(k=0)^∞)((c[k])^2)(((2k+1))^2)

These results are in perfect agreement:

ExpectedEnergy==SpectrallyExpectedEnergy

Had we worked with truncated spectral representations of ψ we would have obtained nearly correct results:

Table[{n,ℰ N[(∑_(k=0)^n)((c[k])^2)(((2k+1))^2),10]},{n,5,50,5}]//TableForm

in which connection we note that the convergence appears to be logarithmically (?) slow.

DIGRESSION: Comparison of direct and spectral derivative structures

Dψ=D[ψ[x],x];
DDψ=D[ψ[x],{x,2}];
DDDψ=D[ψ[x],{x,3}];

Plot[{Dψ,DDψ,DDDψ},{x,0,1}, PlotStyle->{{Red,Thick},{Blue,Thick},{Green,Thick}}]

DSpectralΨ=D[SpectralΨ[x],x];
DDSpectralΨ=D[SpectralΨ[x],{x,2}];
DDDSpectralΨ=D[SpectralΨ[x],{x,3}];

The spectral descriptions of those quite-simple functions are quite complicated, so Mathematica takes a minute or two to construct he following figure:

Fig1=Plot[Evaluate[DSpectralΨ],{x,0,1}, PlotStyle->{Red,Thick}];
Fig2=Plot[Evaluate[DDSpectralΨ],{x,0,1}, PlotStyle->{Blue,Thick}];
Fig3=Plot[Evaluate[DDDSpectralΨ],{x,0,1}, PlotStyle->{Green,Thick}];

Show[Fig1, Fig2, Fig3,PlotRange->{-12,6}]

And the simplest of the functions was evidently posed for Mathematica the greatest computational challenge (notice how the green curve is broken). 

Mathematica has less difficulty if one works in finite approximation; i.e., with truncated series:

Ψ[x_]:=(∑_(k=0)^10)c[k]e[x,2k+1]
DΨ=D[Ψ[x],x];
DDΨ=D[Ψ[x],{x,2}];
DDDΨ=D[Ψ[x],{x,3}];

fig1=Plot[Evaluate[DΨ],{x,0,1}, PlotStyle->{Red}];
fig2=Plot[Evaluate[DDΨ],{x,0,1}, PlotStyle->{Blue}];
fig3=Plot[Evaluate[DDDΨ],{x,0,1}, PlotStyle->{Green}];

Show[fig1, fig2,fig3,PlotRange->{-12,6}]

But, as the figure illustrates, truncation leads to imprecisions that become rapidly more conspicuous as the order of the differentiation increases. I increase the number of terms in the truncated series, and look in finer detail to the results:

ΨΨ[x_]:=(∑_(k=0)^50)c[k]e[x,2k+1]
DΨΨ=D[ΨΨ[x],x];
DDΨΨ=D[ΨΨ[x],{x,2}];
DDDΨΨ=D[ΨΨ[x],{x,3}];

Plot[Evaluate[DΨΨ],{x,0,1}, PlotStyle->{Red}, PlotRange->{-12,6}]

Plot[Evaluate[DDΨΨ],{x,0,1}, PlotStyle->{Blue}, PlotPoints->50, PlotRange->{-13,0}]

Gibbs phenomenon is the dominant feature of the preceding figure. It is even more conspicuous in the following figure

Plot[Evaluate[DDDΨΨ],{x,0,1}, PlotStyle->{Green}, PlotPoints->50, PlotRange->{-60,60}]

of which this is a 10× magnification:

Plot[Evaluate[DDDΨΨ],{x,0,1}, PlotStyle->{Green}, PlotPoints->50, PlotRange->{-600,600}]

If we call that function f(x) then it is only in an approximate (meaning averaged) sense that it is the zero function, though we can expect to have

except possibly for Gibbsian end effects.

The reason one must keep more and more terms to get good approximations to exact results emerges when one compares

((8(Sqrt[30]))/((π^3)(n^3)))e[x,n];
D[%,x]
D[%,x]
D[%,x]
D[%,x]
D[%,x]

The short of it: the (n^(-3)) factors that were responsible for the initially rapid convergence lose force with each differentiation, disappear altogether after 3 differentiations, and thereafter become "divergence factors."

Expected value of energysquared

ExpectedEnergySquared=(∫_0^1)ψ[x]((ℋ^2) D[ψ[x],{x,4}])\[DifferentialD]x

That result is surprising, but understandable: ψ[x] is a quadratic function of x, so is killed by 4th order differentiation. 

If, on the other hand, we proceed from the spectral representation of ψ we already got

SpectrallyExpectedEnergy=ℰ(∑_(k=0)^∞)((c[k])^2)(((2k+1))^2)

and now get

SpectrallyExpectedEnergySquared=(ℰ^2)(∑_(k=0)^∞)((c[k])^2)(((2k+1))^4)

which is non-zero. We have arrived at a contradiction:

ExpectedEnergySquared-SpectrallyExpectedEnergySquared

In third order (as in all higher orders) the results yielded by the two methods differ even more profoundly: the first method yields 0 in every case, while the second method leads in every case to series that (for the reason stated just above) diverge: we have

SpectrallyExpectedEnergyCubed=(ℰ^3)(∑_(k=0)^p)((c[k])^2)(((2k+1))^6)

which blows up in the limit p ⟶ ∞.

The evident source of the contradiction: we have arrived at a situation in which

derivative of infinite sum ≠infinite sum of derivatives

REMARK: It is clear already on dimensional grounds that relaxation of our assumption that the box width a = 1 entails

Variance

By definition

The spectral method supplies

SpectrallyExpectedEnergySquared-(((SpectrallyExpectedEnergy))^2)

while direct calculation leads to the nonsensical result

ExpectedEnergySquared-(((ExpectedEnergy))^2)

Dynamical manifestation of the same problem

Working as before from

H=-ℋ(∂_(x,x))\[Placeholder]

we expect to have

Clear[ℋ]
ϕ[x,t]=(∑_(k=0)^∞)(1/(k!))(((-(I/ℏ)ℋ t))^k)D[ϕ[x],{x,2k}]
ℋ=((ℏ^2)/(2m));

Looking to the motion of the individual energy eigenfunctions, we have

Table[{k,((D[e[x,n],{x,2k}])/(e[x,n]))},{k,0,4}]//TableForm

whence

(∑_(k=0)^∞)(1/(k!))(((-(I/ℏ)ℋ t))^k)(((-(n^2) (π^2)))^k)(nth.eigenstate)

which—when we recall that

ℰ (n^2)

is just (Subscript[ℰ, n])—gives back the familiar fact that the energy eigenstates simply buzz:

(Subscript[(nth.eigenstate), t])=(E^((I/ℏ)(Subscript[ℰ, n])t))(Subscript[(nth.eigenstate), 0])

Here I automate the demonstration:

Clear[ℋ,ℰ]
(∑_(k=0)^∞)(1/(k!))(((-I ℋ t))^k)D[e[x,n],{x,2k}]
%/.(∂_({x,2 k}))((Sqrt[2]) Sin[n π x])->(((-(n^2) (π^2)))^k)(Sqrt[2]) Sin[n π x];
%/.I (n^2) (π^2) t ℋ->(I/ℏ)(Subscript[ℰ, n])t
ℋ=((ℏ^2)/(2m));
ℰ=(n^2) (π^2) ℋ;

Given the spectral representation of a state, we have only to launch each of its components into motion to set the whole state into motion. Which we now demonstrate as it pertains to the parabolic test state. Define

e[x_,n_,t_]:=e[x,n](E^(-I (n^2)t))

where as a matter of computational convenience I have set ℰ/ℏ = 1. Construct the truncated spectral series (I want only provide proof of method, though we know that inclusion of higher-order terms would have little quantitative effect)

Ψ[x_,t_]:=(∑_(k=0)^2)c[k]e[x,2k+1,t]

and use

CompConj[A_] := A /.Complex[0,n_]->-Complex[0,n]

to construct

ConjΨ[x_,t_]:=CompConj[Ψ[x,t]]

Those functions look like this

Ψ[x,t]
ConjΨ[x,t]

and become rapidly more unwieldly as one moves to higher order. Construct

MovingProbDensity[x_,t_]:=ComplexExpand[Ψ[x,t]CompConj[Ψ[x,t]]]

MovingProbDensity[x,t]

TrigReduce[%]

ReducedMovingProbDensity[x_,t_]:=-(1/(759375 (π^6)))64 (-11406979+11390625 Cos[2 π x]+15625 Cos[6 π x]+729 Cos[10 π x]+3375 Cos[16 t-8 π x]+91125 Cos[24 t-6 π x]+421875 Cos[8 t-4 π x]-91125 Cos[24 t-4 π x]-421875 Cos[8 t-2 π x]-3375 Cos[16 t-2 π x]-421875 Cos[8 t+2 π x]-3375 Cos[16 t+2 π x]+421875 Cos[8 t+4 π x]-91125 Cos[24 t+4 π x]+91125 Cos[24 t+6 π x]+3375 Cos[16 t+8 π x])

Note that the first cycle of this periodic function is complete at time t = 2π/8.

I construct now a figure illustrating how the moving density compares with its initial structure

((ψ[x])^2)

Manipulate[Plot[{((ψ[x])^2),Evaluate[ReducedMovingProbDensity[x,t]]},{x,0,1},PlotRange->{0,2.5}],{t,0,((2π)/8)}]

We now abandon the spectral representation, and work directly from

ψ[x]

We then find

Clear[ℋ]
(∑_(k=0)^∞)(1/(k!))(((-I ℋ t))^k)D[ψ[x],{x,2k}]

which by

Table[(∂_({x,2 k}))((Sqrt[30]) (1-x) x),{k,0,10}]

gives

(∑_(k=0)^10)(((((-I t ℋ))^k) (∂_({x,2 k}))((Sqrt[30]) (1-x) x))/(k!))

So we have

ψ[x_,t_]:=(Sqrt[30]) (1-x) x+2 I (Sqrt[30]) t ℋ
Conjψ[x_,t_]:=CompConj[ψ[x,t]]

giving

MMovingProbDensity[x_,t_]:=ComplexExpand[ψ[x,t]CompConj[ψ[x,t]]]

MMovingProbDensity[x,t]

The motion of ψ[x, t] is clearly not periodic, nor is it even unitary:

(∫_0^1)%\[DifferentialD]x

Our Hamiltonian has revcaled itself to be non-Hermitian.

MEMORANDA

I had stewed about this problem for a couple of days when, on January 7th, I brought the problem to the attention of David Griffiths, who promptly responded as follows:

Nick,

Here's my first guess as to your error: H^2 is not a valid operator (or rather, your second wave function is not in the space of functions for which it is well-defined/Hermitian).  Both your wave functions have kinks at each end, which means that the first derivative is discontinuous and the second carries a delta function.  You get away with ignoring this when you calculate <H> for the "accidental" reason reason that it is multiplied by a function that is zero at the ends, so the delta function does not contribute.  But when you go to H^2 you need two more derivatives (hence the second derivative of the delta function).  Integration by parts transfers them over to the other function, which means you now have the product of two deltas.  For some reason that works in the case of the ground state, but not in the case of the the quadratic.

Is this even vaguely relevant to your problem?

David

When (on January 11th—yesterday) I mentioned the problem with Darrell Schroeter he made an essentially similar remark: he suggested that my wave function might better be written

where S[x, a] has unit value inside the box and vanishes outside the box. Those remarks/suggestions have motivated the preceeding discussion. 

I end with record of some commands that figured in earlier versions of this notebook, which I preserve for whatever they may someday prove to be worth:

CoeffList=Table[(∫_0^1)x e[x,n]\[DifferentialD]x,{n,1,20}]

(∑_(k=1)^20)CoeffList[[k]]e[x,k]

Plot[{x,((2 Sin[π x])/π)-((Sin[2 π x])/π)+((2 Sin[3 π x])/(3 π))-((Sin[4 π x])/(2 π))+((2 Sin[5 π x])/(5 π))-((Sin[6 π x])/(3 π))+((2 Sin[7 π x])/(7 π))-((Sin[8 π x])/(4 π))+((2 Sin[9 π x])/(9 π))-((Sin[10 π x])/(5 π))+((2 Sin[11 π x])/(11 π))-((Sin[12 π x])/(6 π))+((2 Sin[13 π x])/(13 π))-((Sin[14 π x])/(7 π))+((2 Sin[15 π x])/(15 π))-((Sin[16 π x])/(8 π))+((2 Sin[17 π x])/(17 π))-((Sin[18 π x])/(9 π))+((2 Sin[19 π x])/(19 π))-((Sin[20 π x])/(10 π))},{x,0,1}]

Plot[-(∑_(k=1)^100)(((-1))^k)(2/(k π))Sin[k π x],{x,0,1}]

Plot[-(∑_(k=1)^100)(((-1))^k)(2/(k π))Sin[k π x],{x,-1.5,1.5}]